
****** programs for bootstrap

*** first stage + reduced form
cap program drop reg_nchild3_edulev3
program define reg_nchild3_edulev3, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild3 fine3b `*' if twinhh == 0
	predict nchild3_hat, xb
	g pscore_d_x = nchild3_hat - _b[fine3b] * fine3b
	drop nchild3_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild3 twinhh fine3b_twinhh fine3b  c.pscore_d_x_dm#c.twinhh c.pscore_d_x_dm#c.fine3b c.pscore_d_x_dm#c.fine3b#c.twinhh pscore_d_x_dm

	scalar alpha_z = _b[twinhh]
	scalar alpha_zx = _b[fine3b_twinhh]
	scalar alpha_x = _b[fine3b]

	scalar r2_fs = e(r2)

	drop pscore_d_x*

* reduced-form
	* pscore_y_x, a summary control for X 
	qui reg edulev3 fine3b `*' if twinhh == 0
	predict edulev3_hat, xb
	g pscore_y_x = edulev3_hat - _b[fine3b] * fine3b
	drop edulev3_hat

	qui sum pscore_y_x, d
	g pscore_y_x_dm = pscore_y_x - r(p50)

	* regression
	reg edulev3 twinhh fine3b_twinhh fine3b c.pscore_y_x_dm#c.twinhh c.pscore_y_x_dm#c.fine3b c.pscore_y_x_dm#c.fine3b#c.twinhh pscore_y_x_dm

	scalar rho_z = _b[twinhh]
	scalar rho_zx = _b[fine3b_twinhh]
	scalar rho_x = _b[fine3b]

	scalar r2_rf = e(r2)

	drop pscore_y_x*
	
* QQ effects
	scalar beta_A = rho_z/alpha_z

	scalar beta_B = rho_zx/alpha_zx

	matrix coef = (alpha_z, alpha_zx, alpha_x, rho_z, rho_zx, rho_x, beta_A, beta_B, r2_fs, r2_rf)
	matrix colnames coef = alpha_z alpha_zx alpha_x rho_z rho_zx rho_x beta_A beta_B r2_fs r2_rf

	ereturn post coef

end


* F-statistic
cap program drop reg_nchild3
program define reg_nchild3, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild3 fine3b `*' if twinhh == 0
	predict nchild3_hat, xb
	g pscore_d_x = nchild3_hat - _b[fine3b] * fine3b
	drop nchild3_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild3 twinhh fine3b_twinhh fine3b  c.pscore_d_x_dm#c.twinhh c.pscore_d_x_dm#c.fine3b c.pscore_d_x_dm#c.fine3b#c.twinhh pscore_d_x_dm

	drop pscore_d_x*
end



*** first stage + reduced form, nchild4
cap program drop reg_nchild4_edulev3
program define reg_nchild4_edulev3, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild4 fine4b `*' if twinhh == 0
	predict nchild4_hat, xb
	g pscore_d_x = nchild4_hat - _b[fine4b] * fine4b
	drop nchild4_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild4 twinhh fine4b_twinhh fine4b  c.pscore_d_x_dm#c.twinhh c.pscore_d_x_dm#c.fine4b c.pscore_d_x_dm#c.fine4b#c.twinhh pscore_d_x_dm

	scalar alpha_z = _b[twinhh]
	scalar alpha_zx = _b[fine4b_twinhh]
	scalar alpha_x = _b[fine4b]

	scalar r2_fs = e(r2)

	drop pscore_d_x*

* reduced-form
	* pscore_y_x, a summary control for X 
	qui reg edulev3 fine4b `*' if twinhh == 0
	predict edulev3_hat, xb
	g pscore_y_x = edulev3_hat - _b[fine4b] * fine4b
	drop edulev3_hat

	qui sum pscore_y_x, d
	g pscore_y_x_dm = pscore_y_x - r(p50)

	* regression
	reg edulev3 twinhh fine4b_twinhh fine4b c.pscore_y_x_dm#c.twinhh c.pscore_y_x_dm#c.fine4b c.pscore_y_x_dm#c.fine4b#c.twinhh pscore_y_x_dm

	scalar rho_z = _b[twinhh]
	scalar rho_zx = _b[fine4b_twinhh]
	scalar rho_x = _b[fine4b]

	scalar r2_rf = e(r2)

	drop pscore_y_x*
	
* QQ effects
	scalar beta_A = rho_z/alpha_z

	scalar beta_B = rho_zx/alpha_zx

	matrix coef = (alpha_z, alpha_zx, alpha_x, rho_z, rho_zx, rho_x, beta_A, beta_B, r2_fs, r2_rf)
	matrix colnames coef = alpha_z alpha_zx alpha_x rho_z rho_zx rho_x beta_A beta_B r2_fs r2_rf

	ereturn post coef

end


* F-statistic
cap program drop reg_nchild4
program define reg_nchild4, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild4 fine4b `*' if twinhh == 0
	predict nchild4_hat, xb
	g pscore_d_x = nchild4_hat - _b[fine4b] * fine4b
	drop nchild4_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild4 twinhh fine4b_twinhh fine4b  c.pscore_d_x_dm#c.twinhh c.pscore_d_x_dm#c.fine4b c.pscore_d_x_dm#c.fine4b#c.twinhh pscore_d_x_dm

	drop pscore_d_x*
end




*** first stage + reduced form, nchild5
cap program drop reg_nchild5_edulev3
program define reg_nchild5_edulev3, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild5 fine5b `*' if twinhh == 0
	predict nchild5_hat, xb
	g pscore_d_x = nchild5_hat - _b[fine5b] * fine5b
	drop nchild5_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild5 twinhh fine5b_twinhh fine5b  c.pscore_d_x_dm#c.twinhh c.pscore_d_x_dm#c.fine5b c.pscore_d_x_dm#c.fine5b#c.twinhh pscore_d_x_dm

	scalar alpha_z = _b[twinhh]
	scalar alpha_zx = _b[fine5b_twinhh]
	scalar alpha_x = _b[fine5b]

	scalar r2_fs = e(r2)

	drop pscore_d_x*

* reduced-form
	* pscore_y_x, a summary control for X 
	qui reg edulev3 fine5b `*' if twinhh == 0
	predict edulev3_hat, xb
	g pscore_y_x = edulev3_hat - _b[fine5b] * fine5b
	drop edulev3_hat

	qui sum pscore_y_x, d
	g pscore_y_x_dm = pscore_y_x - r(p50)

	* regression
	reg edulev3 twinhh fine5b_twinhh fine5b c.pscore_y_x_dm#c.twinhh c.pscore_y_x_dm#c.fine5b c.pscore_y_x_dm#c.fine5b#c.twinhh pscore_y_x_dm

	scalar rho_z = _b[twinhh]
	scalar rho_zx = _b[fine5b_twinhh]
	scalar rho_x = _b[fine5b]

	scalar r2_rf = e(r2)

	drop pscore_y_x*
	
* QQ effects
	scalar beta_A = rho_z/alpha_z

	scalar beta_B = rho_zx/alpha_zx

	matrix coef = (alpha_z, alpha_zx, alpha_x, rho_z, rho_zx, rho_x, beta_A, beta_B, r2_fs, r2_rf)
	matrix colnames coef = alpha_z alpha_zx alpha_x rho_z rho_zx rho_x beta_A beta_B r2_fs r2_rf

	ereturn post coef

end


* F-statistic
cap program drop reg_nchild5
program define reg_nchild5, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild5 fine5b `*' if twinhh == 0
	predict nchild5_hat, xb
	g pscore_d_x = nchild5_hat - _b[fine5b] * fine5b
	drop nchild5_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild5 twinhh fine5b_twinhh fine5b  c.pscore_d_x_dm#c.twinhh c.pscore_d_x_dm#c.fine5b c.pscore_d_x_dm#c.fine5b#c.twinhh pscore_d_x_dm

	drop pscore_d_x*
end






*** first stage + reduced form, fine3b3, 3rd-order polynomial
cap program drop reg_nchild3_edulev3_fine3b3
program define reg_nchild3_edulev3_fine3b3, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild3 fine3b fine3b2 fine3b3 `*' if twinhh == 0
	predict nchild3_hat, xb
	g pscore_d_x = nchild3_hat - _b[fine3b] * fine3b - _b[fine3b2] * fine3b2 - _b[fine3b3] * fine3b3
	drop nchild3_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild3 twinhh fine3b_twinhh fine3b2_twinhh fine3b3_twinhh   fine3b fine3b2 fine3b3 ///
		c.pscore_d_x_dm#c.twinhh  c.pscore_d_x_dm#c.fine3b c.pscore_d_x_dm#c.fine3b2 c.pscore_d_x_dm#c.fine3b3 ///
		c.pscore_d_x_dm#c.fine3b#c.twinhh  c.pscore_d_x_dm#c.fine3b2#c.twinhh c.pscore_d_x_dm#c.fine3b3#c.twinhh  pscore_d_x_dm

	scalar alpha_z = _b[twinhh]

	scalar alpha_zx1 = _b[fine3b_twinhh]
	scalar alpha_zx2 = _b[fine3b2_twinhh]
	scalar alpha_zx3 = _b[fine3b3_twinhh]

	scalar alpha_x1 = _b[fine3b]
	scalar alpha_x2 = _b[fine3b2]
	scalar alpha_x3 = _b[fine3b3]

	scalar r2_fs = e(r2)

	drop pscore_d_x*

* reduced-form
	* pscore_y_x, a summary control for X 
	qui reg edulev3 fine3b fine3b2 fine3b3 `*' if twinhh == 0
	predict edulev3_hat, xb
	g pscore_y_x = edulev3_hat - _b[fine3b] * fine3b - _b[fine3b2] * fine3b2 - _b[fine3b3] * fine3b3
	drop edulev3_hat

	qui sum pscore_y_x, d
	g pscore_y_x_dm = pscore_y_x - r(p50)

	* regression
	reg edulev3 twinhh fine3b_twinhh fine3b2_twinhh fine3b3_twinhh   fine3b fine3b2 fine3b3 ///
		c.pscore_y_x_dm#c.twinhh  c.pscore_y_x_dm#c.fine3b c.pscore_y_x_dm#c.fine3b2 c.pscore_y_x_dm#c.fine3b3 ///
		c.pscore_y_x_dm#c.fine3b#c.twinhh  c.pscore_y_x_dm#c.fine3b2#c.twinhh c.pscore_y_x_dm#c.fine3b3#c.twinhh  pscore_y_x_dm

	scalar rho_z = _b[twinhh]
	
	scalar rho_zx1 = _b[fine3b_twinhh]
	scalar rho_zx2 = _b[fine3b2_twinhh]
	scalar rho_zx3 = _b[fine3b3_twinhh]

	scalar rho_x1 = _b[fine3b]
	scalar rho_x2 = _b[fine3b2]
	scalar rho_x3 = _b[fine3b3]

	scalar r2_rf = e(r2)

	drop pscore_y_x*
	
* QQ effects
	scalar beta_A = rho_z/alpha_z

	scalar rho_zx_20 = rho_zx1 * 0.20 + rho_zx2 * 0.20^2 + rho_zx3 * 0.20^3
	scalar alpha_zx_20 = alpha_zx1 * 0.20 + alpha_zx2 * 0.20^2 + alpha_zx3 * 0.20^3
	scalar beta_B_20 = rho_zx_20 / alpha_zx_20

	scalar rho_zx_30 = rho_zx1 * 0.30 + rho_zx2 * 0.30^2 + rho_zx3 * 0.30^3
	scalar alpha_zx_30 = alpha_zx1 * 0.30 + alpha_zx2 * 0.30^2 + alpha_zx3 * 0.30^3
	scalar beta_B_30 = rho_zx_30 / alpha_zx_30

	scalar rho_zx_40 = rho_zx1 * 0.40 + rho_zx2 * 0.40^2 + rho_zx3 * 0.40^3
	scalar alpha_zx_40 = alpha_zx1 * 0.40 + alpha_zx2 * 0.40^2 + alpha_zx3 * 0.40^3
	scalar beta_B_40 = rho_zx_40 / alpha_zx_40

	scalar rho_zx_50 = rho_zx1 * 0.50 + rho_zx2 * 0.50^2 + rho_zx3 * 0.50^3
	scalar alpha_zx_50 = alpha_zx1 * 0.50 + alpha_zx2 * 0.50^2 + alpha_zx3 * 0.50^3
	scalar beta_B_50 = rho_zx_50 / alpha_zx_50	

	scalar rho_zx_60 = rho_zx1 * 0.60 + rho_zx2 * 0.60^2 + rho_zx3 * 0.60^3
	scalar alpha_zx_60 = alpha_zx1 * 0.60 + alpha_zx2 * 0.60^2 + alpha_zx3 * 0.60^3
	scalar beta_B_60 = rho_zx_60 / alpha_zx_60

	scalar rho_zx_70 = rho_zx1 * 0.70 + rho_zx2 * 0.70^2 + rho_zx3 * 0.70^3
	scalar alpha_zx_70 = alpha_zx1 * 0.70 + alpha_zx2 * 0.70^2 + alpha_zx3 * 0.70^3
	scalar beta_B_70 = rho_zx_70 / alpha_zx_70

	scalar rho_zx_80 = rho_zx1 * 0.80 + rho_zx2 * 0.80^2 + rho_zx3 * 0.80^3
	scalar alpha_zx_80 = alpha_zx1 * 0.80 + alpha_zx2 * 0.80^2 + alpha_zx3 * 0.80^3
	scalar beta_B_80 = rho_zx_80 / alpha_zx_80

	scalar rho_zx_90 = rho_zx1 * 0.90 + rho_zx2 * 0.90^2 + rho_zx3 * 0.90^3
	scalar alpha_zx_90 = alpha_zx1 * 0.90 + alpha_zx2 * 0.90^2 + alpha_zx3 * 0.90^3
	scalar beta_B_90 = rho_zx_90 / alpha_zx_90

	scalar rho_zx_100 = rho_zx1 * 1.00 + rho_zx2 * 1.00^2 + rho_zx3 * 1.00^3
	scalar alpha_zx_100 = alpha_zx1 * 1.00 + alpha_zx2 * 1.00^2 + alpha_zx3 * 1.00^3
	scalar beta_B_100 = rho_zx_100 / alpha_zx_100

	scalar rho_zx_110 = rho_zx1 * 1.10 + rho_zx2 * 1.10^2 + rho_zx3 * 1.10^3
	scalar alpha_zx_110 = alpha_zx1 * 1.10 + alpha_zx2 * 1.10^2 + alpha_zx3 * 1.10^3
	scalar beta_B_110 = rho_zx_110 / alpha_zx_110

	scalar rho_zx_120 = rho_zx1 * 1.20 + rho_zx2 * 1.20^2 + rho_zx3 * 1.20^3
	scalar alpha_zx_120 = alpha_zx1 * 1.20 + alpha_zx2 * 1.20^2 + alpha_zx3 * 1.20^3
	scalar beta_B_120 = rho_zx_120 / alpha_zx_120

	matrix coef = (alpha_z, alpha_zx1, alpha_zx2, alpha_zx3, alpha_x1, alpha_x2, alpha_x3, ///
		rho_z, rho_zx1, rho_zx2, rho_zx3, rho_x1, rho_x2, rho_x3, beta_A, r2_fs, r2_rf, ///
		beta_B_20, beta_B_30, beta_B_40, beta_B_50, beta_B_60, beta_B_70, ///
		beta_B_80, beta_B_90, beta_B_100, beta_B_110, beta_B_120)
	matrix colnames coef = alpha_z alpha_zx1 alpha_zx2 alpha_zx3 alpha_x1 alpha_x2 alpha_x3 ///
		rho_z rho_zx1 rho_zx2 rho_zx3 rho_x1 rho_x2 rho_x3 beta_A r2_fs r2_rf ///
		beta_B_20 beta_B_30 beta_B_40 beta_B_50 beta_B_60 beta_B_70 ///
		beta_B_80 beta_B_90 beta_B_100 beta_B_110 beta_B_120

	ereturn post coef

end






*** first stage + reduced form, fine3b4, 4th-order polynomial
cap program drop reg_nchild3_edulev3_fine3b4
program define reg_nchild3_edulev3_fine3b4, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild3 fine3b fine3b2 fine3b3 fine3b4 `*' if twinhh == 0
	predict nchild3_hat, xb
	g pscore_d_x = nchild3_hat - _b[fine3b] * fine3b - _b[fine3b2] * fine3b2 - _b[fine3b3] * fine3b3 - _b[fine3b4] * fine3b4
	drop nchild3_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild3 twinhh fine3b_twinhh fine3b2_twinhh fine3b3_twinhh fine3b4_twinhh  fine3b fine3b2 fine3b3 fine3b4 ///
		c.pscore_d_x_dm#c.twinhh  c.pscore_d_x_dm#c.fine3b c.pscore_d_x_dm#c.fine3b2 c.pscore_d_x_dm#c.fine3b3 c.pscore_d_x_dm#c.fine3b4 ///
		c.pscore_d_x_dm#c.fine3b#c.twinhh  c.pscore_d_x_dm#c.fine3b2#c.twinhh c.pscore_d_x_dm#c.fine3b3#c.twinhh c.pscore_d_x_dm#c.fine3b4#c.twinhh  pscore_d_x_dm

	scalar alpha_z = _b[twinhh]

	scalar alpha_zx1 = _b[fine3b_twinhh]
	scalar alpha_zx2 = _b[fine3b2_twinhh]
	scalar alpha_zx3 = _b[fine3b3_twinhh]
	scalar alpha_zx4 = _b[fine3b4_twinhh]

	scalar alpha_x1 = _b[fine3b]
	scalar alpha_x2 = _b[fine3b2]
	scalar alpha_x3 = _b[fine3b3]
	scalar alpha_x4 = _b[fine3b4]

	scalar r2_fs = e(r2)

	drop pscore_d_x*

* reduced-form
	* pscore_y_x, a summary control for X 
	qui reg edulev3 fine3b fine3b2 fine3b3 fine3b4 `*' if twinhh == 0
	predict edulev3_hat, xb
	g pscore_y_x = edulev3_hat - _b[fine3b] * fine3b - _b[fine3b2] * fine3b2 - _b[fine3b3] * fine3b3 - _b[fine3b4] * fine3b4
	drop edulev3_hat

	qui sum pscore_y_x, d
	g pscore_y_x_dm = pscore_y_x - r(p50)

	* regression
	reg edulev3 twinhh fine3b_twinhh fine3b2_twinhh fine3b3_twinhh fine3b4_twinhh  fine3b fine3b2 fine3b3 fine3b4 ///
		c.pscore_y_x_dm#c.twinhh  c.pscore_y_x_dm#c.fine3b c.pscore_y_x_dm#c.fine3b2 c.pscore_y_x_dm#c.fine3b3 c.pscore_y_x_dm#c.fine3b4 ///
		c.pscore_y_x_dm#c.fine3b#c.twinhh  c.pscore_y_x_dm#c.fine3b2#c.twinhh c.pscore_y_x_dm#c.fine3b3#c.twinhh c.pscore_y_x_dm#c.fine3b4#c.twinhh  pscore_y_x_dm

	scalar rho_z = _b[twinhh]
	
	scalar rho_zx1 = _b[fine3b_twinhh]
	scalar rho_zx2 = _b[fine3b2_twinhh]
	scalar rho_zx3 = _b[fine3b3_twinhh]
	scalar rho_zx4 = _b[fine3b4_twinhh]

	scalar rho_x1 = _b[fine3b]
	scalar rho_x2 = _b[fine3b2]
	scalar rho_x3 = _b[fine3b3]
	scalar rho_x4 = _b[fine3b4]

	scalar r2_rf = e(r2)

	drop pscore_y_x*
	
* QQ effects
	scalar beta_A = rho_z/alpha_z

	scalar rho_zx_20 = rho_zx1 * 0.20 + rho_zx2 * 0.20^2 + rho_zx3 * 0.20^3 + rho_zx4 * 0.20^4
	scalar alpha_zx_20 = alpha_zx1 * 0.20 + alpha_zx2 * 0.20^2 + alpha_zx3 * 0.20^3 + alpha_zx4 * 0.20^4
	scalar beta_B_20 = rho_zx_20 / alpha_zx_20

	scalar rho_zx_30 = rho_zx1 * 0.30 + rho_zx2 * 0.30^2 + rho_zx3 * 0.30^3 + rho_zx4 * 0.30^4
	scalar alpha_zx_30 = alpha_zx1 * 0.30 + alpha_zx2 * 0.30^2 + alpha_zx3 * 0.30^3 + alpha_zx4 * 0.30^4
	scalar beta_B_30 = rho_zx_30 / alpha_zx_30

	scalar rho_zx_40 = rho_zx1 * 0.40 + rho_zx2 * 0.40^2 + rho_zx3 * 0.40^3 + rho_zx4 * 0.40^4
	scalar alpha_zx_40 = alpha_zx1 * 0.40 + alpha_zx2 * 0.40^2 + alpha_zx3 * 0.40^3 + alpha_zx4 * 0.40^4
	scalar beta_B_40 = rho_zx_40 / alpha_zx_40

	scalar rho_zx_50 = rho_zx1 * 0.50 + rho_zx2 * 0.50^2 + rho_zx3 * 0.50^3 + rho_zx4 * 0.50^4
	scalar alpha_zx_50 = alpha_zx1 * 0.50 + alpha_zx2 * 0.50^2 + alpha_zx3 * 0.50^3 + alpha_zx4 * 0.50^4
	scalar beta_B_50 = rho_zx_50 / alpha_zx_50

	scalar rho_zx_60 = rho_zx1 * 0.60 + rho_zx2 * 0.60^2 + rho_zx3 * 0.60^3 + rho_zx4 * 0.60^4
	scalar alpha_zx_60 = alpha_zx1 * 0.60 + alpha_zx2 * 0.60^2 + alpha_zx3 * 0.60^3 + alpha_zx4 * 0.60^4
	scalar beta_B_60 = rho_zx_60 / alpha_zx_60

	scalar rho_zx_70 = rho_zx1 * 0.70 + rho_zx2 * 0.70^2 + rho_zx3 * 0.70^3 + rho_zx4 * 0.70^4
	scalar alpha_zx_70 = alpha_zx1 * 0.70 + alpha_zx2 * 0.70^2 + alpha_zx3 * 0.70^3 + alpha_zx4 * 0.70^4
	scalar beta_B_70 = rho_zx_70 / alpha_zx_70

	scalar rho_zx_80 = rho_zx1 * 0.80 + rho_zx2 * 0.80^2 + rho_zx3 * 0.80^3 + rho_zx4 * 0.80^4
	scalar alpha_zx_80 = alpha_zx1 * 0.80 + alpha_zx2 * 0.80^2 + alpha_zx3 * 0.80^3 + alpha_zx4 * 0.80^4
	scalar beta_B_80 = rho_zx_80 / alpha_zx_80

	scalar rho_zx_90 = rho_zx1 * 0.90 + rho_zx2 * 0.90^2 + rho_zx3 * 0.90^3 + rho_zx4 * 0.90^4
	scalar alpha_zx_90 = alpha_zx1 * 0.90 + alpha_zx2 * 0.90^2 + alpha_zx3 * 0.90^3 + alpha_zx4 * 0.90^4
	scalar beta_B_90 = rho_zx_90 / alpha_zx_90

	scalar rho_zx_100 = rho_zx1 * 1.00 + rho_zx2 * 1.00^2 + rho_zx3 * 1.00^3 + rho_zx4 * 1.00^4
	scalar alpha_zx_100 = alpha_zx1 * 1.00 + alpha_zx2 * 1.00^2 + alpha_zx3 * 1.00^3 + alpha_zx4 * 1.00^4
	scalar beta_B_100 = rho_zx_100 / alpha_zx_100

	scalar rho_zx_110 = rho_zx1 * 1.10 + rho_zx2 * 1.10^2 + rho_zx3 * 1.10^3 + rho_zx4 * 1.10^4
	scalar alpha_zx_110 = alpha_zx1 * 1.10 + alpha_zx2 * 1.10^2 + alpha_zx3 * 1.10^3 + alpha_zx4 * 1.10^4
	scalar beta_B_110 = rho_zx_110 / alpha_zx_110

	scalar rho_zx_120 = rho_zx1 * 1.20 + rho_zx2 * 1.20^2 + rho_zx3 * 1.20^3 + rho_zx4 * 1.20^4
	scalar alpha_zx_120 = alpha_zx1 * 1.20 + alpha_zx2 * 1.20^2 + alpha_zx3 * 1.20^3 + alpha_zx4 * 1.20^4
	scalar beta_B_120 = rho_zx_120 / alpha_zx_120

	matrix coef = (alpha_z, alpha_zx1, alpha_zx2, alpha_zx3, alpha_zx4, alpha_x1, alpha_x2, alpha_x3, alpha_x4, ///
		rho_z, rho_zx1, rho_zx2, rho_zx3, rho_zx4, rho_x1, rho_x2, rho_x3, rho_x4, beta_A, r2_fs, r2_rf, ///
		beta_B_20, beta_B_30, beta_B_40, beta_B_50, beta_B_60, beta_B_70, ///
		beta_B_80, beta_B_90, beta_B_100, beta_B_110, beta_B_120)
	matrix colnames coef = alpha_z alpha_zx1 alpha_zx2 alpha_zx3 alpha_zx4 alpha_x1 alpha_x2 alpha_x3 alpha_x4 ///
		rho_z rho_zx1 rho_zx2 rho_zx3 rho_zx4 rho_x1 rho_x2 rho_x3 rho_x4 beta_A r2_fs r2_rf ///
		beta_B_20 beta_B_30 beta_B_40 beta_B_50 beta_B_60 beta_B_70 ///
		beta_B_80 beta_B_90 beta_B_100 beta_B_110 beta_B_120

	ereturn post coef

end



*** first stage + reduced form, fine3b5, 5th-order polynomial
cap program drop reg_nchild3_edulev3_fine3b5
program define reg_nchild3_edulev3_fine3b5, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild3 fine3b fine3b2 fine3b3 fine3b4 fine3b5 `*' if twinhh == 0
	predict nchild3_hat, xb
	g pscore_d_x = nchild3_hat - _b[fine3b] * fine3b - _b[fine3b2] * fine3b2 - _b[fine3b3] * fine3b3 - _b[fine3b4] * fine3b4 - _b[fine3b5] * fine3b5
	drop nchild3_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild3 twinhh fine3b_twinhh fine3b2_twinhh fine3b3_twinhh fine3b4_twinhh fine3b5_twinhh  fine3b fine3b2 fine3b3 fine3b4 fine3b5 ///
		c.pscore_d_x_dm#c.twinhh  c.pscore_d_x_dm#c.fine3b c.pscore_d_x_dm#c.fine3b2 c.pscore_d_x_dm#c.fine3b3 c.pscore_d_x_dm#c.fine3b4 c.pscore_d_x_dm#c.fine3b5 ///
		c.pscore_d_x_dm#c.fine3b#c.twinhh  c.pscore_d_x_dm#c.fine3b2#c.twinhh c.pscore_d_x_dm#c.fine3b3#c.twinhh c.pscore_d_x_dm#c.fine3b4#c.twinhh c.pscore_d_x_dm#c.fine3b5#c.twinhh  pscore_d_x_dm

	scalar alpha_z = _b[twinhh]

	scalar alpha_zx1 = _b[fine3b_twinhh]
	scalar alpha_zx2 = _b[fine3b2_twinhh]
	scalar alpha_zx3 = _b[fine3b3_twinhh]
	scalar alpha_zx4 = _b[fine3b4_twinhh]
	scalar alpha_zx5 = _b[fine3b5_twinhh]

	scalar alpha_x1 = _b[fine3b]
	scalar alpha_x2 = _b[fine3b2]
	scalar alpha_x3 = _b[fine3b3]
	scalar alpha_x4 = _b[fine3b4]
	scalar alpha_x5 = _b[fine3b5]

	scalar r2_fs = e(r2)

	drop pscore_d_x*

* reduced-form
	* pscore_y_x, a summary control for X 
	qui reg edulev3 fine3b fine3b2 fine3b3 fine3b4 fine3b5 `*' if twinhh == 0
	predict edulev3_hat, xb
	g pscore_y_x = edulev3_hat - _b[fine3b] * fine3b - _b[fine3b2] * fine3b2 - _b[fine3b3] * fine3b3 - _b[fine3b4] * fine3b4 - _b[fine3b5] * fine3b5
	drop edulev3_hat

	qui sum pscore_y_x, d
	g pscore_y_x_dm = pscore_y_x - r(p50)

	* regression
	*reg edulev3 twinhh fine3b_twinhh fine3b2_twinhh fine3b3_twinhh fine3b4_twinhh fine3b5_twinhh ///
	*	fine3b fine3b2 fine3b3 fine3b4 fine3b5  c.pscore_y_x_dm#c.twinhh pscore_y_x_dm
	reg edulev3 twinhh fine3b_twinhh fine3b2_twinhh fine3b3_twinhh fine3b4_twinhh fine3b5_twinhh  fine3b fine3b2 fine3b3 fine3b4 fine3b5 ///
		c.pscore_y_x_dm#c.twinhh  c.pscore_y_x_dm#c.fine3b c.pscore_y_x_dm#c.fine3b2 c.pscore_y_x_dm#c.fine3b3 c.pscore_y_x_dm#c.fine3b4 c.pscore_y_x_dm#c.fine3b5 ///
		c.pscore_y_x_dm#c.fine3b#c.twinhh  c.pscore_y_x_dm#c.fine3b2#c.twinhh c.pscore_y_x_dm#c.fine3b3#c.twinhh c.pscore_y_x_dm#c.fine3b4#c.twinhh c.pscore_y_x_dm#c.fine3b5#c.twinhh  pscore_y_x_dm

	scalar rho_z = _b[twinhh]
	
	scalar rho_zx1 = _b[fine3b_twinhh]
	scalar rho_zx2 = _b[fine3b2_twinhh]
	scalar rho_zx3 = _b[fine3b3_twinhh]
	scalar rho_zx4 = _b[fine3b4_twinhh]
	scalar rho_zx5 = _b[fine3b5_twinhh]

	scalar rho_x1 = _b[fine3b]
	scalar rho_x2 = _b[fine3b2]
	scalar rho_x3 = _b[fine3b3]
	scalar rho_x4 = _b[fine3b4]
	scalar rho_x5 = _b[fine3b5]

	scalar r2_rf = e(r2)

	drop pscore_y_x*
	
* QQ effects
	scalar beta_A = rho_z/alpha_z

	scalar rho_zx_20 = rho_zx1 * 0.20 + rho_zx2 * 0.20^2 + rho_zx3 * 0.20^3 + rho_zx4 * 0.20^4 + rho_zx5 * 0.20^5
	scalar alpha_zx_20 = alpha_zx1 * 0.20 + alpha_zx2 * 0.20^2 + alpha_zx3 * 0.20^3 + alpha_zx4 * 0.20^4 + alpha_zx5 * 0.20^5
	scalar beta_B_20 = rho_zx_20 / alpha_zx_20

	scalar rho_zx_30 = rho_zx1 * 0.30 + rho_zx2 * 0.30^2 + rho_zx3 * 0.30^3 + rho_zx4 * 0.30^4 + rho_zx5 * 0.30^5
	scalar alpha_zx_30 = alpha_zx1 * 0.30 + alpha_zx2 * 0.30^2 + alpha_zx3 * 0.30^3 + alpha_zx4 * 0.30^4 + alpha_zx5 * 0.30^5
	scalar beta_B_30 = rho_zx_30 / alpha_zx_30

	scalar rho_zx_40 = rho_zx1 * 0.40 + rho_zx2 * 0.40^2 + rho_zx3 * 0.40^3 + rho_zx4 * 0.40^4 + rho_zx5 * 0.40^5
	scalar alpha_zx_40 = alpha_zx1 * 0.40 + alpha_zx2 * 0.40^2 + alpha_zx3 * 0.40^3 + alpha_zx4 * 0.40^4 + alpha_zx5 * 0.40^5
	scalar beta_B_40 = rho_zx_40 / alpha_zx_40

	scalar rho_zx_50 = rho_zx1 * 0.50 + rho_zx2 * 0.50^2 + rho_zx3 * 0.50^3 + rho_zx4 * 0.50^4 + rho_zx5 * 0.50^5
	scalar alpha_zx_50 = alpha_zx1 * 0.50 + alpha_zx2 * 0.50^2 + alpha_zx3 * 0.50^3 + alpha_zx4 * 0.50^4 + alpha_zx5 * 0.50^5
	scalar beta_B_50 = rho_zx_50 / alpha_zx_50

	scalar rho_zx_60 = rho_zx1 * 0.60 + rho_zx2 * 0.60^2 + rho_zx3 * 0.60^3 + rho_zx4 * 0.60^4 + rho_zx5 * 0.60^5
	scalar alpha_zx_60 = alpha_zx1 * 0.60 + alpha_zx2 * 0.60^2 + alpha_zx3 * 0.60^3 + alpha_zx4 * 0.60^4 + alpha_zx5 * 0.60^5
	scalar beta_B_60 = rho_zx_60 / alpha_zx_60

	scalar rho_zx_70 = rho_zx1 * 0.70 + rho_zx2 * 0.70^2 + rho_zx3 * 0.70^3 + rho_zx4 * 0.70^4 + rho_zx5 * 0.70^5
	scalar alpha_zx_70 = alpha_zx1 * 0.70 + alpha_zx2 * 0.70^2 + alpha_zx3 * 0.70^3 + alpha_zx4 * 0.70^4 + alpha_zx5 * 0.70^5
	scalar beta_B_70 = rho_zx_70 / alpha_zx_70

	scalar rho_zx_80 = rho_zx1 * 0.80 + rho_zx2 * 0.80^2 + rho_zx3 * 0.80^3 + rho_zx4 * 0.80^4 + rho_zx5 * 0.80^5
	scalar alpha_zx_80 = alpha_zx1 * 0.80 + alpha_zx2 * 0.80^2 + alpha_zx3 * 0.80^3 + alpha_zx4 * 0.80^4 + alpha_zx5 * 0.80^5
	scalar beta_B_80 = rho_zx_80 / alpha_zx_80

	scalar rho_zx_90 = rho_zx1 * 0.90 + rho_zx2 * 0.90^2 + rho_zx3 * 0.90^3 + rho_zx4 * 0.90^4 + rho_zx5 * 0.90^5
	scalar alpha_zx_90 = alpha_zx1 * 0.90 + alpha_zx2 * 0.90^2 + alpha_zx3 * 0.90^3 + alpha_zx4 * 0.90^4 + alpha_zx5 * 0.90^5
	scalar beta_B_90 = rho_zx_90 / alpha_zx_90

	scalar rho_zx_100 = rho_zx1 * 1.00 + rho_zx2 * 1.00^2 + rho_zx3 * 1.00^3 + rho_zx4 * 1.00^4 + rho_zx5 * 1.00^5
	scalar alpha_zx_100 = alpha_zx1 * 1.00 + alpha_zx2 * 1.00^2 + alpha_zx3 * 1.00^3 + alpha_zx4 * 1.00^4 + alpha_zx5 * 1.00^5
	scalar beta_B_100 = rho_zx_100 / alpha_zx_100

	scalar rho_zx_110 = rho_zx1 * 1.10 + rho_zx2 * 1.10^2 + rho_zx3 * 1.10^3 + rho_zx4 * 1.10^4 + rho_zx5 * 1.10^5
	scalar alpha_zx_110 = alpha_zx1 * 1.10 + alpha_zx2 * 1.10^2 + alpha_zx3 * 1.10^3 + alpha_zx4 * 1.10^4 + alpha_zx5 * 1.10^5
	scalar beta_B_110 = rho_zx_110 / alpha_zx_110

	scalar rho_zx_120 = rho_zx1 * 1.20 + rho_zx2 * 1.20^2 + rho_zx3 * 1.20^3 + rho_zx4 * 1.20^4 + rho_zx5 * 1.20^5
	scalar alpha_zx_120 = alpha_zx1 * 1.20 + alpha_zx2 * 1.20^2 + alpha_zx3 * 1.20^3 + alpha_zx4 * 1.20^4 + alpha_zx5 * 1.20^5
	scalar beta_B_120 = rho_zx_120 / alpha_zx_120


	matrix coef = (alpha_z, alpha_zx1, alpha_zx2, alpha_zx3, alpha_zx4, alpha_zx5, alpha_x1, alpha_x2, alpha_x3, alpha_x4, alpha_x5, ///
		rho_z, rho_zx1, rho_zx2, rho_zx3, rho_zx4, rho_zx5, rho_x1, rho_x2, rho_x3, rho_x4, rho_x5, beta_A, r2_fs, r2_rf, ///
		beta_B_20, beta_B_30, beta_B_40, beta_B_50, beta_B_60, beta_B_70, ///
		beta_B_80, beta_B_90, beta_B_100, beta_B_110, beta_B_120)
	matrix colnames coef = alpha_z alpha_zx1 alpha_zx2 alpha_zx3 alpha_zx4 alpha_zx5 alpha_x1 alpha_x2 alpha_x3 alpha_x4 alpha_x5 ///
		rho_z rho_zx1 rho_zx2 rho_zx3 rho_zx4 rho_zx5 rho_x1 rho_x2 rho_x3 rho_x4 rho_x5 beta_A r2_fs r2_rf ///
		beta_B_20 beta_B_30 beta_B_40 beta_B_50 beta_B_60 beta_B_70 ///
		beta_B_80 beta_B_90 beta_B_100 beta_B_110 beta_B_120

	ereturn post coef

end



*** first stage + reduced form, male as the shifter
cap program drop reg_nchild3_edulev3_bymale
program define reg_nchild3_edulev3_bymale, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild3 male `*' if twinhh == 0
	predict nchild3_hat, xb
	g pscore_d_x = nchild3_hat - _b[male] * male
	drop nchild3_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild3 twinhh male_twinhh male  c.pscore_d_x_dm#c.twinhh  c.pscore_d_x_dm#c.male c.pscore_d_x_dm#c.male#c.twinhh pscore_d_x_dm

	scalar alpha_z = _b[twinhh]
	scalar alpha_zx = _b[male_twinhh]
	scalar alpha_x = _b[male]

	scalar r2_fs = e(r2)

	drop pscore_d_x*

* reduced-form
	* pscore_y_x, a summary control for X 
	qui reg edulev3 male `*' if twinhh == 0
	predict edulev3_hat, xb
	g pscore_y_x = edulev3_hat - _b[male] * male
	drop edulev3_hat

	qui sum pscore_y_x, d
	g pscore_y_x_dm = pscore_y_x - r(p50)

	* regression
	reg edulev3 twinhh male_twinhh male c.pscore_y_x_dm#c.twinhh c.pscore_y_x_dm#c.male c.pscore_y_x_dm#c.male#c.twinhh  pscore_y_x_dm

	scalar rho_z = _b[twinhh]
	scalar rho_zx = _b[male_twinhh]
	scalar rho_x = _b[male]

	scalar r2_rf = e(r2)

	drop pscore_y_x*
	
* QQ effects
	scalar beta_A = rho_z/alpha_z

	scalar beta_B = rho_zx/alpha_zx

	matrix coef = (alpha_z, alpha_zx, alpha_x, rho_z, rho_zx, rho_x, beta_A, beta_B, r2_fs, r2_rf)
	matrix colnames coef = alpha_z alpha_zx alpha_x rho_z rho_zx rho_x beta_A beta_B r2_fs r2_rf

	ereturn post coef

end


* F-statistic
cap program drop reg_nchild3_bymale
program define reg_nchild3_bymale, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild3 male `*' if twinhh == 0
	predict nchild3_hat, xb
	g pscore_d_x = nchild3_hat - _b[male] * male
	drop nchild3_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild3 twinhh male_twinhh male  c.pscore_d_x_dm#c.twinhh  c.pscore_d_x_dm#c.male c.pscore_d_x_dm#c.male#c.twinhh pscore_d_x_dm

	drop pscore_d_x*
end



*** first stage + reduced form, eduyP2dm as the shifter
cap program drop reg_nchild3_edulev3_byeduyP2dm
program define reg_nchild3_edulev3_byeduyP2dm, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild3 eduyP2dm `*' if twinhh == 0
	predict nchild3_hat, xb
	g pscore_d_x = nchild3_hat - _b[eduyP2dm] * eduyP2dm
	drop nchild3_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild3 twinhh eduyP2dm_twinhh eduyP2dm  c.pscore_d_x_dm#c.twinhh c.pscore_d_x_dm#c.eduyP2dm c.pscore_d_x_dm#c.eduyP2dm#c.twinhh  pscore_d_x_dm

	scalar alpha_z = _b[twinhh]
	scalar alpha_zx = _b[eduyP2dm_twinhh]
	scalar alpha_x = _b[eduyP2dm]

	scalar r2_fs = e(r2)

	drop pscore_d_x*

* reduced-form
	* pscore_y_x, a summary control for X 
	qui reg edulev3 eduyP2dm `*' if twinhh == 0
	predict edulev3_hat, xb
	g pscore_y_x = edulev3_hat - _b[eduyP2dm] * eduyP2dm
	drop edulev3_hat

	qui sum pscore_y_x, d
	g pscore_y_x_dm = pscore_y_x - r(p50)

	* regression
	reg edulev3 twinhh eduyP2dm_twinhh eduyP2dm c.pscore_y_x_dm#c.twinhh c.pscore_y_x_dm#c.eduyP2dm c.pscore_y_x_dm#c.eduyP2dm#c.twinhh  pscore_y_x_dm

	scalar rho_z = _b[twinhh]
	scalar rho_zx = _b[eduyP2dm_twinhh]
	scalar rho_x = _b[eduyP2dm]

	scalar r2_rf = e(r2)

	drop pscore_y_x*
	
* QQ effects
	scalar beta_A = rho_z/alpha_z
	*scalar beta_dx = (rho_z + rho_zx) / (alpha_z + alpha_zx) - rho_z/alpha_z

	scalar beta_B = rho_zx/alpha_zx

	matrix coef = (alpha_z, alpha_zx, alpha_x, rho_z, rho_zx, rho_x, beta_A, beta_B, r2_fs, r2_rf)
	matrix colnames coef = alpha_z alpha_zx alpha_x rho_z rho_zx rho_x beta_A beta_B r2_fs r2_rf

	ereturn post coef

end


* F-statistic
cap program drop reg_nchild3_byeduyP2dm
program define reg_nchild3_byeduyP2dm, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild3 eduyP2dm `*' if twinhh == 0
	predict nchild3_hat, xb
	g pscore_d_x = nchild3_hat - _b[eduyP2dm] * eduyP2dm
	drop nchild3_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild3 twinhh eduyP2dm_twinhh eduyP2dm  c.pscore_d_x_dm#c.twinhh c.pscore_d_x_dm#c.eduyP2dm c.pscore_d_x_dm#c.eduyP2dm#c.twinhh  pscore_d_x_dm

	drop pscore_d_x*
end







*** first stage + reduced form, gender-interactive model
cap program drop reg_nchild3_edulev3_gender
program define reg_nchild3_edulev3_gender, eclass

* first-stage
	* pscore_d_x, a summary control for X 
	qui reg nchild3 fine3b `*' if twinhh == 0
	predict nchild3_hat, xb
	g pscore_d_x = nchild3_hat - _b[fine3b] * fine3b
	drop nchild3_hat

	qui sum pscore_d_x, d
	g pscore_d_x_dm = pscore_d_x - r(p50)

	* regression
	reg nchild3 twinhh c.twinhh#c.female fine3b_twinhh c.fine3b_twinhh#c.female fine3b c.fine3b#c.female female ///
	  	c.pscore_d_x_dm#c.twinhh c.pscore_d_x_dm#c.fine3b c.pscore_d_x_dm#c.fine3b#c.twinhh pscore_d_x_dm ///
	  	c.pscore_d_x_dm#c.twinhh#c.female c.pscore_d_x_dm#c.fine3b#c.female c.pscore_d_x_dm#c.fine3b#c.twinhh#c.female c.pscore_d_x_dm#c.female

	scalar alpha_z = _b[twinhh]
	scalar alpha_zx = _b[fine3b_twinhh]
	scalar alpha_x = _b[fine3b]

	scalar alpha_z_f = _b[c.twinhh#c.female]
	scalar alpha_zx_f = _b[c.fine3b_twinhh#c.female]
	scalar alpha_x_f = _b[c.fine3b#c.female]

	scalar alpha_f = _b[female]

	scalar r2_fs = e(r2)

	drop pscore_d_x*

* reduced-form
	* pscore_y_x, a summary control for X 
	qui reg edulev3 fine3b `*' if twinhh == 0
	predict edulev3_hat, xb
	g pscore_y_x = edulev3_hat - _b[fine3b] * fine3b
	drop edulev3_hat

	qui sum pscore_y_x, d
	g pscore_y_x_dm = pscore_y_x - r(p50)

	* regression
	reg edulev3 twinhh c.twinhh#c.female fine3b_twinhh c.fine3b_twinhh#c.female fine3b c.fine3b#c.female female ///
	  	c.pscore_y_x_dm#c.twinhh c.pscore_y_x_dm#c.fine3b c.pscore_y_x_dm#c.fine3b#c.twinhh pscore_y_x_dm ///
	  	c.pscore_y_x_dm#c.twinhh#c.female c.pscore_y_x_dm#c.fine3b#c.female c.pscore_y_x_dm#c.fine3b#c.twinhh#c.female c.pscore_y_x_dm#c.female

	scalar rho_z = _b[twinhh]
	scalar rho_zx = _b[fine3b_twinhh]
	scalar rho_x = _b[fine3b]

	scalar rho_z_f = _b[c.twinhh#c.female]
	scalar rho_zx_f = _b[c.fine3b_twinhh#c.female]
	scalar rho_x_f = _b[c.fine3b#c.female]

	scalar rho_f = _b[female]

	scalar r2_rf = e(r2)

	drop pscore_y_x*
	
* QQ effects
	scalar beta_A = rho_z/alpha_z

	scalar beta_B = rho_zx/alpha_zx

	scalar beta_A_f = (rho_z + rho_z_f) / (alpha_z + alpha_z_f)

	scalar beta_B_f = (rho_zx + rho_zx_f) / (alpha_zx + alpha_zx_f)

	matrix coef = (alpha_z, alpha_zx, alpha_x, alpha_z_f, alpha_zx_f, alpha_x_f, alpha_f, ///
					rho_z, rho_zx, rho_x, rho_z_f, rho_zx_f, rho_x_f, rho_f, ///
					beta_A, beta_B, beta_A_f, beta_B_f, r2_fs, r2_rf)
	matrix colnames coef = alpha_z alpha_zx alpha_x alpha_z_f alpha_zx_f alpha_x_f alpha_f ///
					rho_z rho_zx rho_x rho_z_f rho_zx_f rho_x_f rho_f ///
					beta_A beta_B beta_A_f beta_B_f r2_fs r2_rf

	ereturn post coef

end



